Expoloring data


Filtering counts

After viewing the size factors, it appears that we have 5 samples with really low values that should be removed:

  • 19355
  • 19372
  • 19382
  • 19496
  • 19505


Without filtering low reads, we had a total of 38828 counts. After filtering out the low counts (those with a base mean less than 3), we now have 24384 counts remaining for all C. virginica samples.



PCA plots of Global Gene Expression

# going to use the filter_counts_out object created above because it is filtered for low reads and has outliers removed

## create the DESeq object
dds <- DESeqDataSetFromMatrix(countData = filter_counts_out, 
                              colData = expDesign,
                              design = ~ treat)
dds <- DESeq(dds) # differential expression analysis on gamma-poisson distribution
vsd <-varianceStabilizingTransformation(dds, blind = TRUE) # quickly estimate dispersion trend and apply a variance stabilizing transformation

## Save dds and vsd data 
save(dds, vsd, file = "Data/transformed_counts.RData")
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 1500
## 
## adonis2(formula = pcaData[6:31] ~ pcaData$infect * pcaData$pCO2, permutations = 1500, method = "eu")
##                             Df SumOfSqs      R2      F Pr(>F)
## pcaData$infect               1   1118.3 0.04736 1.1893 0.1799
## pcaData$pCO2                 1    945.3 0.04003 1.0053 0.3817
## pcaData$infect:pcaData$pCO2  1    864.0 0.03659 0.9188 0.5623
## Residual                    22  20686.8 0.87603              
## Total                       25  23614.3 1.00000




Differential gene expression


DEGs

## log2 fold change (MLE): treat S 400 vs N 2800 
## Wald test p-value: treat S 400 vs N 2800 
## DataFrame with 24384 rows and 6 columns
##                        baseMean      log2FoldChange             lfcSE
##                       <numeric>           <numeric>         <numeric>
## LOC111126949   3.45225877157332  -0.942224682074193 0.817783710717942
## LOC111110729  0.485466216780156  -0.512448068281346  1.46364594751556
## LOC111120752   1.69362126099025  -0.282706315518707  0.84681189235072
## LOC111113860   5.93338041130341 -0.0448270626067275 0.601266321422372
## LOC111109550  0.331392203994242  -0.986053856261322  2.48365300058619
## ...                         ...                 ...               ...
## LOC111117112  0.646302530637881    1.67887461151877  1.80337428132596
## LOC111117113   0.14363842589226   0.768137359775011  3.83879218552636
## LOC111117117 0.0778207961311771   0.287240683720725  4.24705363394974
## LOC111116908  0.109745657738521 -0.0817811312489327  4.24579857131686
## LOC111117715  0.621026843813011  -0.445504180182541   1.6799677252913
##                            stat            pvalue              padj
##                       <numeric>         <numeric>         <numeric>
## LOC111126949  -1.15216856208471 0.249251813595227 0.999651667882862
## LOC111110729  -0.35011750563802 0.726250513749595 0.999651667882862
## LOC111120752 -0.333847833352841 0.738494386321806 0.999651667882862
## LOC111113860 -0.074554421243291 0.940569239720562 0.999651667882862
## LOC111109550 -0.397017560838247 0.691354510914985 0.999651667882862
## ...                         ...               ...               ...
## LOC111117112  0.930962933709109 0.351872738132339 0.999651667882862
## LOC111117113   0.20009870882596 0.841403383153489 0.999651667882862
## LOC111117117 0.0676329306097302 0.946077840571377 0.999651667882862
## LOC111116908 -0.019261660645283 0.984632368623089 0.999651667882862
## LOC111117715 -0.265186154159772 0.790866060568095 0.999651667882862
## 
## out of 24380 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up)       : 2, 0.0082%
## LFC < 0 (down)     : 0, 0%
## outliers [1]       : 36, 0.15%
## low counts [2]     : 3, 0.012%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
summary(res_con_sponge)
## 
## out of 24380 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 0, 0%
## LFC < 0 (down)     : 0, 0%
## outliers [1]       : 36, 0.15%
## low counts [2]     : 3, 0.012%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
summary(res_con_pco2)
## 
## out of 24380 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 0, 0%
## LFC < 0 (down)     : 0, 0%
## outliers [1]       : 36, 0.15%
## low counts [2]     : 3, 0.012%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
summary(res_trt)
## 
## out of 24380 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 6, 0.025%
## LFC < 0 (down)     : 5, 0.021%
## outliers [1]       : 36, 0.15%
## low counts [2]     : 11339, 47%
## (mean count < 2)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
#summary(res_trt_pco2)



Venn of DEGs

It is just a demo, not very exciting :D

Session Information

Session information from the last full knit of Rmarkdown on 29 June 2022.

## R version 3.6.3 (2020-02-29)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Catalina 10.15.7
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
##  [1] grid      parallel  stats4    stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] ggvenn_0.1.9                adegenet_2.1.4             
##  [3] ade4_1.7-17                 ggrepel_0.9.1              
##  [5] pdftools_3.0.1              ggpubr_0.4.0               
##  [7] data.table_1.14.0           vegan_2.5-7                
##  [9] lattice_0.20-44             permute_0.9-5              
## [11] DESeq2_1.26.0               SummarizedExperiment_1.16.1
## [13] DelayedArray_0.12.3         BiocParallel_1.20.1        
## [15] matrixStats_0.57.0          Biobase_2.46.0             
## [17] GenomicRanges_1.38.0        GenomeInfoDb_1.22.1        
## [19] IRanges_2.20.2              S4Vectors_0.24.4           
## [21] BiocGenerics_0.32.0         plotly_4.9.4.1             
## [23] arrayQualityMetrics_3.42.0  forcats_0.5.1              
## [25] stringr_1.4.0               dplyr_1.0.7                
## [27] purrr_0.3.4                 readr_2.0.1                
## [29] tidyr_1.1.3                 tibble_3.1.3               
## [31] ggplot2_3.3.5               tidyverse_1.3.1            
## [33] knitr_1.33                 
## 
## loaded via a namespace (and not attached):
##   [1] utf8_1.2.2             tidyselect_1.1.1       RSQLite_2.2.7         
##   [4] AnnotationDbi_1.48.0   htmlwidgets_1.5.3      beadarray_2.36.1      
##   [7] munsell_0.5.0          units_0.7-2            codetools_0.2-18      
##  [10] preprocessCore_1.48.0  withr_2.4.2            colorspace_2.0-2      
##  [13] highr_0.9              rstudioapi_0.13        setRNG_2013.9-1       
##  [16] ggsignif_0.6.2         labeling_0.4.2         GenomeInfoDbData_1.2.2
##  [19] hwriter_1.3.2          farver_2.1.0           bit64_4.0.5           
##  [22] LearnBayes_2.15.1      coda_0.19-4            vctrs_0.3.8           
##  [25] generics_0.1.0         xfun_0.30              qpdf_1.1              
##  [28] R6_2.5.1               illuminaio_0.28.0      gridSVG_1.7-2         
##  [31] locfit_1.5-9.4         bitops_1.0-7           cachem_1.0.5          
##  [34] assertthat_0.2.1       promises_1.2.0.1       scales_1.1.1          
##  [37] nnet_7.3-16            gtable_0.3.0           affy_1.64.0           
##  [40] rlang_0.4.11           genefilter_1.68.0      systemfonts_1.0.2     
##  [43] splines_3.6.3          rstatix_0.7.0          lazyeval_0.2.2        
##  [46] hexbin_1.28.2          broom_0.7.9            checkmate_2.0.0       
##  [49] BiocManager_1.30.16    yaml_2.3.5             reshape2_1.4.4        
##  [52] abind_1.4-5            modelr_0.1.8           crosstalk_1.1.1       
##  [55] backports_1.2.1        httpuv_1.6.2           Hmisc_4.5-0           
##  [58] tools_3.6.3            spData_0.3.10          affyio_1.56.0         
##  [61] ellipsis_0.3.2         raster_3.4-13          jquerylib_0.1.4       
##  [64] RColorBrewer_1.1-2     proxy_0.4-26           Rcpp_1.0.7            
##  [67] plyr_1.8.6             base64enc_0.1-3        zlibbioc_1.32.0       
##  [70] classInt_0.4-3         RCurl_1.98-1.4         deldir_0.2-10         
##  [73] rpart_4.1-15           openssl_1.4.4          haven_2.4.3           
##  [76] cluster_2.1.2          fs_1.5.0               magrittr_2.0.1        
##  [79] openxlsx_4.2.4         gmodels_2.18.1         reprex_2.0.1          
##  [82] hms_1.1.0              mime_0.11              evaluate_0.14         
##  [85] xtable_1.8-4           XML_3.99-0.3           rio_0.5.27            
##  [88] jpeg_0.1-9             gcrma_2.58.0           readxl_1.3.1          
##  [91] gridExtra_2.3          compiler_3.6.3         KernSmooth_2.23-20    
##  [94] crayon_1.4.1           htmltools_0.5.2        spdep_1.1-8           
##  [97] mgcv_1.8-36            later_1.3.0            tzdb_0.1.2            
## [100] Formula_1.2-4          geneplotter_1.64.0     expm_0.999-6          
## [103] lubridate_1.7.10       DBI_1.1.1              dbplyr_2.1.1          
## [106] MASS_7.3-54            boot_1.3-28            sf_1.0-5              
## [109] Matrix_1.3-4           car_3.0-11             cli_3.0.1             
## [112] vsn_3.54.0             gdata_2.18.0           igraph_1.2.6          
## [115] pkgconfig_2.0.3        sp_1.4-5               foreign_0.8-75        
## [118] xml2_1.3.2             svglite_2.0.0          annotate_1.64.0       
## [121] bslib_0.3.1            BeadDataPackR_1.38.0   affyPLM_1.62.0        
## [124] XVector_0.26.0         rvest_1.0.1            digest_0.6.27         
## [127] Biostrings_2.54.0      rmarkdown_2.13         base64_2.0            
## [130] cellranger_1.1.0       htmlTable_2.2.1        curl_4.3.2            
## [133] gtools_3.9.2           shiny_1.7.1            lifecycle_1.0.0       
## [136] nlme_3.1-151           jsonlite_1.7.2         carData_3.0-4         
## [139] seqinr_4.2-8           viridisLite_0.4.0      askpass_1.1           
## [142] limma_3.42.2           fansi_0.5.0            pillar_1.6.2          
## [145] fastmap_1.1.0          httr_1.4.2             survival_3.2-12       
## [148] glue_1.4.2             zip_2.2.0              png_0.1-7             
## [151] bit_4.0.4              class_7.3-19           stringi_1.7.3         
## [154] sass_0.4.0             blob_1.2.2             latticeExtra_0.6-29   
## [157] memoise_2.0.0          e1071_1.7-8            ape_5.5